rm(list=ls())
setwd("/Volumes/Disk2/Future_PM/Dataverse")
library(fields); library(maps); library(ncdf4);library(abind)
source('Function/get_geo.R')
source('Function/get_met.R')
source('Function/functions.R')

ss=load("Data/ACCMIP/mask_high_resolution.Rdata")
lon.in=lon.out
lat.in=lat.out

ss=load("Data/AQS_obs/AQS_species_slopes.Rdata")
mask=sp.dissolve(land.mask,lon.in,lat.in, lon.frame, lat.frame)<0.3
mask[is.na(mask)]=FALSE

cr=0.1
total_slope[total_sig> cr]=NA;total_slope[mask]=NA
SO4_slope[SO4_sig> cr]=NA; SO4_slope[mask]=NA
NH4_slope[NH4_sig> cr]=NA; NH4_slope[mask]=NA
NO3_slope[NO3_sig> cr]=NA; NO3_slope[mask]=NA
OC_slope[OC_sig> cr]=NA; OC_slope[mask]=NA
total_slope1=total_slope
SO4_slope1= SO4_slope
NH4_slope1= NH4_slope
NO3_slope1= NO3_slope
OC_slope1= OC_slope

#ss2=load("Data/ACCMIP/GEOS-Chem_species_slopes.Rdata")
ss2=load("Data/GEOS-Chem/GEOS-Chem_species_slopes.Rdata")
mask=sp.dissolve(land.mask,lon.in,lat.in, lon.out, lat.out)<0.3
mask[is.na(mask)]=FALSE

total_slope[total_sig> cr]=NA;total_slope[mask]=NA
SO4_slope[SO4_sig> cr]=NA; SO4_slope[mask]=NA
NH4_slope[NH4_sig> cr]=NA; NH4_slope[mask]=NA
NO3_slope[NO3_sig> cr]=NA; NO3_slope[mask]=NA
OC_slope[OC_sig> cr]=NA; OC_slope[mask]=NA
total_slope2=total_slope
SO4_slope2= SO4_slope
NH4_slope2= NH4_slope
NO3_slope2= NO3_slope
OC_slope2= OC_slope

mai=c(0.1, 0.25, 0.2, 0.05)
mgp=c(1.2, 0.4, 0)
tcl=-0.2
ps=12
lon.ind=3:27

dev.new(width=5,height=6)
close.screen(all.screens = TRUE)
m <- rbind(c(0.00, 0.47, 0.8, 1), c(0.43, 0.90, 0.8, 1), c(0.61, 0.98, 0.8, 1),
		c(0.00, 0.47, 0.62, 0.82),c(0.42, 0.89, 0.62, 0.82),c(0.60, 0.97, 0.62, 0.82),
		c(0.00, 0.47, 0.44, 0.64),c(0.42, 0.89, 0.44, 0.64),c(0.60, 0.97, 0.44, 0.64),
		c(0.00, 0.47, 0.26, 0.46),c(0.42, 0.89, 0.26, 0.46),c(0.60, 0.97, 0.26, 0.46),
  		 c(0.00, 0.47, 0.08, 0.28),c(0.42, 0.89, 0.08, 0.28),c(0.60, 0.97, 0.08, 0.28)					
				  )
				  
split.screen(m)
#----- total slope --------------------------
screen(1)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.frame[lon.ind],lat.frame, total_slope1[lon.ind,],zlim=c(-2.5,2.5),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
mtext("dPM2.5/dT",side=2,font=1,cex=1.00)
title("AQS",cex.main=1.0, font.main=1)
text("(a)",x=-123,y=27)
box()

screen(2)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.out[lon.ind],lat.out, total_slope2[lon.ind,],zlim=c(-2.5,2.5),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
box()
title("GEOS-Chem",cex.main=1.0,font.main=1)
text("(b)",x=-123,y=28)

screen(3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image.plot(legend.shrink=0.9,legend.only=TRUE,zlim=c(-2.5,2.5),col=rwb.colors(32),horizontal=FALSE,legend.width=2.0,axis.args = list(mgp = c(0, 0.3, 0),tcl=-0.2,padj=0.45))

#----- SO4 --------------------------
screen(4)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.frame[lon.ind],lat.frame, SO4_slope1[lon.ind,],zlim=c(-1,1),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
box()
mtext("dSO4/dT",side=2,font=1,cex=1.00)
text("(c)",x=-123,y=27)

screen(5)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.out[lon.ind],lat.out, SO4_slope2[lon.ind,],zlim=c(-1,1),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
box()
text("(d)",x=-123,y=28)

screen(6)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image.plot(legend.shrink=0.9,legend.only=TRUE,zlim=c(-1,1),col=rwb.colors(32),horizontal=FALSE,legend.width=2.0,axis.args = list(mgp = c(0, 0.3, 0),tcl=-0.2,padj=0.45))

#----- NH4 --------------------------
screen(7+3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.frame[lon.ind],lat.frame, NH4_slope1[lon.ind,],zlim=c(-1,1),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
box()
mtext("dNH4/dT",side=2,font=1,cex=1.00)
text("(g)",x=-123,y=27)

screen(8+3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.out[lon.ind],lat.out, NH4_slope2[lon.ind,],zlim=c(-1,1),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
box()
text("(h)",x=-123,y=28)

screen(9+3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image.plot(legend.shrink=0.9,legend.only=TRUE,zlim=c(-1,1),col=rwb.colors(32),horizontal=FALSE,legend.width=2.0,axis.args = list(mgp = c(0, 0.3, 0),tcl=-0.2,padj=0.45))


#----- NO3 --------------------------
screen(10+3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.frame[lon.ind],lat.frame, NO3_slope1[lon.ind,],zlim=c(-1,1),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
box()
mtext("dNO3/dT",side=2,font=1,cex=1.00)
text("(i)",x=-123,y=27)

screen(11+3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.out[lon.ind],lat.out, NO3_slope2[lon.ind,],zlim=c(-1,1),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
box()
text("(j)",x=-123,y=28)

screen(12+3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image.plot(NO3_slope,legend.shrink=0.9,legend.only=TRUE,zlim=c(-1,1),col=rwb.colors(32),horizontal=FALSE,legend.width=2.0,axis.args = list(mgp = c(0, 0.3, 0),tcl=-0.2,padj=0.45))


#----- OA --------------------------
screen(13-6)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.frame[lon.ind],lat.frame, 1.8*OC_slope1[lon.ind,],zlim=c(-1,1),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
box()
mtext("dOA/dT",side=2,font=1,cex=1.00)
text("(e)",x=-123,y=27)

screen(14-6)
OC_slope2[OC_slope2>=1]=1; OC_slope2[OC_slope2<=-1]=-1
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.out[lon.ind],lat.out, OC_slope2[lon.ind,],zlim=c(-1,1),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
box()
text("(f)",x=-123,y=28)

screen(15-6)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image.plot(legend.shrink=0.9,legend.only=TRUE,zlim=c(-1,1),col=rwb.colors(32),horizontal=FALSE,legend.width=2.0,axis.args = list(mgp = c(0, 0.3, 0),tcl=-0.2,padj=0.45))